      program gadget_writer
      
      implicit none
      include 'mpif.h'
      integer*8 int4max
      parameter (int4max=2147483647)
 
      integer npart
      integer nmax
      parameter (nmax=20000000)

C     Input file
      real*8 xx(3,nmax)
      real dxx(3,nmax),ms(nmax)
      real*8 lbox
      integer*4 ilist(nmax)
      common /abc/ xx,dxx,ms,ilist
      integer nmfiles, nwhich
      integer*8 first_index,itotal
      integer*4 ifirst
      real*8  redshift_input,omega0_input,lambda0_input
      real omega,omegaz,lambda,zinit,zics,ugas
      real pcentre(3),transcale,rsphcut,sph_cut,split_factor,tgas,baryon_frac
      integer ishift,itrunc,iendian,isph
      integer num_parallel_files,ncount
      integer nijk,nrat
      integer iunit
      integer ngas_total,ngas_local,nextra
      integer nsteps,istep
      integer id8flag
      integer*4 nhighword(6)
      integer*4 nlowword(6),dum
      real*4 mspmin(6)
      real*4 mspmax(6)

C     Output GADGET file

      integer*4 np(6),nall(6)
      integer*8 nall8(6)                
      real*8 massarr(6)     
      real*8 expansion,redshift 
      integer*4 flagsfr,flagfeedback,flagcooling
      integer*4 NumFiles
      real*8 BoxSize
      real*8 Omega0,OmegaLambda,HubbleParam         
      character unused(256-6*4-6*4-6*8-2*8-4*4-4*8) 
      real aexp1_temp,aexp2_temp,lingro
      integer*4 inum1,inum2

      integer*4 ntall(6)
      integer*8 ntall8(6)

C     
      real rhocrit
      parameter(rhocrit=27.755)
      real mcut(6),mmcut(6),ssph_cut,gas_shift

C     Timing information

      real vfac,dfac,gfactor1,gfactor2
      integer nspecies
      
      double precision mtotal
      real opz
      double precision facc
      parameter(facc=1.0e-5)

      integer i,j,k
      integer nlres
            
      character*512 paramfile
      character*512 infilebase,infile,outfilebase,outfile

      integer iargc
      external iargc

      logical fexist

      integer IERR, MYID, NUMPROCS
      integer MYIDMINUS, MYIDPLUS
      integer status(MPI_STATUS_SIZE)
      integer stag,rtag,NUMID

      real*4 xx0,yy0,zz0,rr0

      call MPI_INIT(IERR)

      call MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS,IERR)

      call MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)


      if (MYID.eq.0) then !----------------------------------------

      if(iargc().eq.1) then
         call getarg(1,paramfile)

         open (12,file=paramfile,status='old')
         read (12,*) infilebase
         read (12,*) outfilebase
         read (12,*) ishift 
         read (12,*) pcentre(1),pcentre(2),pcentre(3)
         read (12,*) itrunc
         read (12,*) rsphcut 
         read (12,*) zinit
         read (12,*) iendian
         read (12,*) isph
         read (12,*) baryon_frac  
         read (12,*) sph_cut
         read (12,*) split_factor
         read (12,*) tgas  
         read (12,*) nspecies
         print*, 'Segregating particles into ',nspecies,' species'
         do i=1,nspecies-1
         read (12,*) mcut(i)
         if (mcut(i).ge.0) then
             print*,'Mass cut must be negative! '
             stop
         endif
         enddo
         if (ishift.eq.0) then 
              transcale = 0.0
              print*,'Not translating particle positions'
         else
              transcale = 1.0
              print*,'Translating particle positions by ',(pcentre(j),j=1,3)
         endif
         
         if ((iendian.lt.0).or.(iendian.gt.2)) then
            print*,'Value of iendian out of range 0-2'
            stop
         endif

         if ((isph.lt.0).or.(isph.gt.1)) then
            print*,'Value of isph out of range 0-1'
            stop
         endif

         if ((split_factor.lt.0.0).or.(split_factor.gt.1.0)) then
            print*,'Value of split_factor  out of range 0-1 (inc)'
            stop
         endif



         close(12)

         call flush(6)
     
         call file_name_generate(infilebase,MYID,infile)

         print*,'Reading file ',infile
         call flush(6)

         open (10,file=infile,form='unformatted',status='old')
         read (10) npart,BoxSize,num_parallel_files
         close(10) 
         

      else
         print*, 'Usage: format_converter.x paramfile'
         goto 666
      end if
      

      endif !------------- end processor 0 -----------------------------

c  Communicate parameters to all processors:

      call MPI_BCAST(infilebase,132,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(outfilebase,132,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(ishift,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(pcentre,3,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(itrunc,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(rsphcut,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(zinit,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(iendian,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(isph,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(sph_cut,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(baryon_frac,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(split_factor,1,MPI_REAL,0,MPI_COMM_WORLD,IERR) 
      call MPI_BCAST(tgas,1,MPI_REAL,0,MPI_COMM_WORLD,IERR) 
      call MPI_BCAST(nspecies,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(mcut,6,MPI_REAL,0,MPI_COMM_WORLD,IERR)   
      call MPI_BCAST(transcale,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)  
      call MPI_BCAST(num_parallel_files,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)  

      do j=1,6
        ntall(j) = 0  ! Accumulate number of particles of each species
        mspmin(j) = 1.e30
        mspmax(j) = -1.e30
      enddo 

      nsteps = num_parallel_files/NUMPROCS 
      if (nsteps*NUMPROCS.ne.num_parallel_files) nsteps=nsteps+1

c   Loop over files to convert
      
      do istep = 1,nsteps
      if (MYID.eq.0) print*,'Starting step ', istep,'/',nsteps

      nijk = NUMPROCS * (istep-1) +  MYID

      if (nijk.lt. num_parallel_files) then        

      call file_name_generate(infilebase,nijk,infile)

         open (10,file=infile,form='unformatted',status='unknown')

         print*, 'Reading ',infile
         call flush(6)

         read (10) npart,lbox,nmfiles,nwhich,first_index,
     &     redshift_input,omega0_input,lambda0_input,itotal
         read (10) xx0,yy0,zz0,rr0

         print*,'MYID,nwhich,nijk',MYID,nwhich,nijk

         if (itotal.gt.int4max) then
                  id8flag = 1 
         else
                  id8flag = 0
         endif


         if(npart.gt.nmax) stop 'Need to increase Nmax'
         print*, 'Read header...  MYID=',MYID
         print*, 'Reading ',npart,' particles MYID = ',MYID
         call flush(6)

         omega = omega0_input
         lambda= lambda0_input
         zics  = redshift_input
         if (nmfiles.ne.num_parallel_files) stop 'nmfiles.ne.num_parallel_files'

         read (10)  (xx(1,i),i=1,npart)
         read (10)  (xx(2,i),i=1,npart)
         read (10)  (xx(3,i),i=1,npart)
         write(*,*) 'Read positions...'

         read (10) (dxx(1,i),i=1,npart)
         read (10) (dxx(2,i),i=1,npart)
         read (10) (dxx(3,i),i=1,npart)
         write(*,*) 'Read displacements...'

         read (10) (ms(i),i=1,npart)
         write(*,*) 'Read masses...'
         write(*,*)minval(ms(1:npart)),maxval(ms(1:npart)),sum(ms(1:npart))

         close (10)
         call flush(6)

         mtotal=0.0
 
         print*,'Here!'
         print*,'nspecies = ',nspecies
         call flush(6)
         do i=1,nspecies-1
         print*,i,omega,rhocrit,lbox
         call flush(6)
         print*,'mcut(i) = ',mcut(i)
         call flush(6)
           mmcut(i)=(10.0**mcut(i))*omega*rhocrit*(lbox**3.0)
         enddo 
         print*,'Here2!'

         ssph_cut = (10.0**sph_cut)*omega*rhocrit*(lbox**3.0)
         
         mmcut(nspecies) = omega*rhocrit*(lbox**3.0)
 
         print*,'Initial displacements prepared for redshift ',zics
         print*,'Value of zinit read in from parameter file ',zinit

c   Recale displacements using linear growth factors-----------------------

         aexp1_temp  = 1./(1.+zics)
         aexp2_temp  = 1./(1.+zinit)
         print*,'omega,lambda ',omega,lambda
         print*,'aexp1_temp,aexp2_temp',aexp1_temp,aexp2_temp
         gfactor1    = lingro(aexp1_temp,omega,lambda)
         gfactor2    = lingro(aexp2_temp,omega,lambda)
         print*,'gfactor1,gfactor2 ',gfactor1,gfactor2
         dfac = gfactor2/gfactor1


c OLD VERSION          dfac = (1.+zics)/(1.+zinit)
         print*,'Scaling displacements by factor: ',dfac 
         print*,'WARNING: assuming Omega very close to 1!!!'
c---------------------------------------------------------------------------

c  Uses approximation dlnD/dlna = Omega**(6/11) which is the asymptotic
c  limit as Omega tends to one for flat universe with a cosmological 
c  constant. This factor is accurate to better than 0.1% for Omega>0.7, 
c  and in error by 1.5% for Omega=0.25.
c  For an open universe the same limit gives Omega**(4/7)  
      
         opz=1.0+zinit             ! Scaling factor
         omegaz=omega*(opz**3.0)/(omega*(opz**3.0)+lambda+
     & (1.0-omega-lambda)*opz**2.0)
         vfac=100.0*(omegaz**0.54545454)*sqrt(omega*(opz**3.0)+lambda+
     & (1.0-omega-lambda)*opz**2.0)/sqrt(opz)
         
         write(*,*) vfac

         do i=1,npart
            do j=1,3
               xx(j,i)=lbox*xx(j,i)
               dxx(j,i)=dfac*lbox*dxx(j,i)
            end do
         end do

         do i=1,npart

            do j=1,3
               xx(j,i)=xx(j,i)+dxx(j,i)              ! Compute positions
               dxx(j,i)=vfac*dxx(j,i)                ! Transform to velocity

c   Periodic boundary conditions apply
               if (xx(j,i).lt.0.) xx(j,i)=xx(j,i)+lbox
               if (xx(j,i).ge.lbox) xx(j,i)=xx(j,i)-lbox
               

            end do

c  Include coordinate translation if needed 

            do j=1,3
               xx(j,i)=xx(j,i)-transcale*pcentre(j)   ! Compute positions
               if (xx(j,i).lt.0.) xx(j,i)=xx(j,i)+lbox
               if (xx(j,i).ge.lbox) xx(j,i)=xx(j,i)-lbox
            end do




            mtotal=mtotal+ms(i)                      ! Compute total mass
            ms(i)=ms(i)*omega*rhocrit*lbox**3.0      ! Scale mass

         end do


c     Now cut-out unwanted particles

      if (itrunc.ne.0) then
      j=0
      do i=1,npart
       if ((xx(1,i)**2+xx(2,i)**2+xx(3,i)**2).le.rsphcut*rsphcut) then
        j=j+1
        xx(1,j) = xx(1,i)
        xx(2,j) = xx(2,i)
        xx(3,j) = xx(3,i)
        dxx(1,j)=dxx(1,i)
        dxx(2,j)=dxx(2,i)
        dxx(3,j)=dxx(3,i)
        ms(j) = ms(i)
       endif
      enddo
      
      npart = j
      print*,'New value of npart = ',npart
      endif



         if(abs(mtotal-1.0).gt.facc) then
            write(*,*) 'Warning! Mass inaccuracy'
            write(*,*) 'Mass ratio:',mtotal
         end if

      do i=1,6
         np(i)=0
         nall(i)=0
         massarr(i)=0.0
      end do

C    Set to zero - values should be overwritten by the parameter file.

      flagsfr=0
      flagcooling=0
      flagfeedback=0
C
      NumFiles=num_parallel_files



c   Find number of particles of different species - no assumed ordering
c   
      ngas_local = 0.0

      do j=1,npart

      if ((isph.eq.1).and.(ms(j).lt.ssph_cut)) ngas_local=ngas_local+1
      if (nspecies.eq.1) then
            np(2) = np(2)+1
            mspmin(2) = min(mspmin(2),ms(j))
            mspmax(2) = max(mspmax(2),ms(j))
      else 
      do i=1,nspecies-1
         if (ms(j).le.mmcut(i)) then
            np(i+1) = np(i+1)+1
            mspmin(i+1) = min(mspmin(i+1),ms(j))
            mspmax(i+1) = max(mspmax(i+1),ms(j))
            goto 41
         endif
      enddo
 41   continue
      if (ms(j).gt.mmcut(nspecies-1)) then
            np(nspecies+1)=np(nspecies+1)+1
            mspmin(nspecies+1) = min(mspmin(nspecies+1),ms(j))
            mspmax(nspecies+1) = max(mspmax(nspecies+1),ms(j))
      endif                 
      endif

      enddo


      print*,'MYID,ngas_local',MYID,ngas_local
      call flush(6)

      np(1) = ngas_local




      do j=1,6
        ntall(j) = ntall(j) + np(j)
      enddo
      
      do j=1,6
         print*,'MYID,species: ', MYID, j,np(j),ntall(j)
      end do


      print*,'MYID, first_index = ',MYID,first_index
  

      if (isph.eq.1) then
c------------------------- Add gas particles ---------------------------      
c First make room
      do i=npart,1,-1
       do j=1,3
         xx(j,i+ngas_local) =  xx(j,i)
        dxx(j,i+ngas_local) = dxx(j,i)
       enddo
         ms(i+ngas_local)   =  ms(i)
      enddo

c   Refine npart
      npart = npart + ngas_local

      print*,'lbox = ',lbox



c   Now loop over dark matter particles and spawn gas particles.
      ncount = 0
      do i=ngas_local+1,npart
        if (ms(i).lt.ssph_cut) then
           ncount = ncount + 1
           gas_shift = 0.5*(ms(i)/(omega*rhocrit))**(1./3.)
           do k=1,3
             xx(k,ncount) =  xx(k,i) + (1.-baryon_frac)*gas_shift*split_factor
            dxx(k,ncount) = dxx(k,i)
             xx(k,i)      =  xx(k,i) -  baryon_frac*gas_shift*split_factor
           enddo

             ms(ncount) = baryon_frac*ms(i)
             mspmin(1) = min(mspmin(1),ms(ncount))
             mspmax(1) = max(mspmax(1),ms(ncount))
             ms(i)      = (1.0-baryon_frac)*ms(i)
        endif
      enddo

      if (ncount.ne.ngas_local) stop 'Error spawning gas particles'

      do i=1,npart  ! Make sure particles lie within 0<=x,y,z<lbox
       do k=1,3
         if (xx(k,i).lt.0.0)  xx(k,i) = xx(k,i) + lbox
         if (xx(k,i).ge.lbox) xx(k,i) = xx(k,i) - lbox
       enddo
      enddo
          

c-----------------------------------------------------------------------
      endif
      
      redshift=zinit
      Omega0=omega
      OmegaLambda=lambda
      BoxSize=lbox
      HubbleParam=0.719

c    Sort by mass

      call indexx(npart,ms,ilist)

      print*,'Sorted particles ',nijk


      expansion=1.0/(1.0+redshift)

      write(*,'(a18,F6.4)') 'Expansion factor: ',expansion
      write(*,'(a8,F6.4,a15,F6.4)') 'Omega0 :',Omega0,
     &                      ' OmegaLambda :',OmegaLambda
      write(*,'(a18,F6.4)') 'Hubble Parameter : ',HubbleParam


      call file_name_generate(outfilebase,nwhich,outfile)
      write(*,'(a10,a50)')  'Writing: ',outfile

      call flush(6)
      print*,MYID,nwhich,outfile
      do j=1,6
       print*,'Species ',j,'number,mmin,mmax',np(j),mspmin(j),mspmax(j)
      enddo
      call flush(6)
   

      
      iunit = iendian + 1  ! Choose unit number - Makefile used to select endianness
                           ! of the output data.
      ugas = 0.0           ! Set gas temperature to zero.

      open(iunit,file=outfile,status='unknown',form='unformatted')

      write(6,*)' np = ',np
      write(6,*)' nall = ',nall
      write(6,*)' ntall = ',ntall
      call flush(6)

      write(iunit) np,massarr,expansion,redshift,flagsfr,flagfeedback
     &       ,nall,flagcooling,NumFiles,BoxSize,Omega0,OmegaLambda
     &       ,HubbleParam,unused 

      write(*,*) 'Writing positions...'
      write(iunit) ((real(xx(j,ilist(i))),j=1,3),i=1,npart)
      write(*,*) 'Writing velocities...'      
      write(iunit) ((dxx(j,ilist(i)),j=1,3),i=1,npart)
      write(*,*) 'Writing identities...'  
      if (id8flag.eq.0) then
        ifirst = first_index    
        write(iunit) (i+ifirst,i=0,npart-1)        ! Integer*4 indexes
      else
        write(iunit) (i+first_index,i=0,npart-1)  ! Integer*8 indexes
      endif
      write(*,*) 'Writing masses...'
      write(iunit) (ms(ilist(i)),i=1,npart)
      if (isph.eq.1) write (iunit) (ugas,i=1,ngas_local) 

      close(iunit)
      call flush(6)

      endif   ! If processor is assigned a file to work on.
      enddo ! loop over steps


      do j=1,6
         ntall8(j) = ntall(j)   ! Transfer to integer*8 variables
      enddo

      call MPI_ALLreduce(ntall8,nall8,6,MPI_INTEGER8,MPI_SUM,MPI_COMM_WORLD,IERR)

      if (MYID.eq.0) then
        do j=1,6
            nhighword(j) = int(nall8(j)/(int4max+1))
            nlowword(j)  = int(nall8(j)-(int4max+1)*nhighword(j))
            print*,'nall(',j,') = ', nall8(j)
            print*,'j,nhigh,nlow ',j,nhighword(j),nlowword(j)
        enddo

      print*,'Now writing in nall values '
      do j=0,num_parallel_files-1
        call file_name_generate(outfilebase,j,outfile)
        open (iunit,file=outfile,access='direct',recl=4,form='unformatted')
        write(6,*)' file = ',trim(outfile)
        write(6,*)' nlowword  = ',nlowword
        write(6,*)' nhighword = ',nhighword
        call flush(6)

        write (iunit,rec=26) nlowword(1)
        write (iunit,rec=27) nlowword(2)
        write (iunit,rec=28) nlowword(3)
        write (iunit,rec=29) nlowword(4)
        write (iunit,rec=30) nlowword(5)
        write (iunit,rec=31) nlowword(6)

        write (iunit,rec=44) nhighword(1)
        write (iunit,rec=45) nhighword(2)
        write (iunit,rec=46) nhighword(3)
        write (iunit,rec=47) nhighword(4)
        write (iunit,rec=48) nhighword(5)
        write (iunit,rec=49) nhighword(6)

        !NB skip block 50, this is flag_entropy_instead_u
        write (iunit,rec=51) real(xx0)
        write (iunit,rec=52) real(yy0)
        write (iunit,rec=53) real(zz0)
        write (iunit,rec=54) real(rr0)

        close (iunit)
       enddo
        print*,'Written nall to all files '
        call flush(6)
      endif



 666  continue

      call MPI_FINALIZE(ierr)

      end



c-------------------- Linear growth rates  - by V.R.Eke -------------

      function lingro(a,omega,lambda)
c
c     To calculate the linear growth factor D(a) for different cosmological 
c     models. Normalised such that D(1) = 1. (Doesn't include closed models 
c     or lambda models where omega+lambda isn't one.)
c
      implicit none
      real a,omega,lambda,func0,func1,func2,x,int,aofx,aofxn
      real xn,w,dn,lingro,sum
      external func0,func1,func2

      w = omega**(-1.0) - 1.0
      sum = omega + lambda
      if (sum.gt.1 .or. omega.le.0 .or.sum.ne.1.and.lambda.gt.0) then
         write(*,*) 'Cannot cope with this cosmology!'
         print*,sum,omega,lambda
         stop
      endif
      if (omega .eq. 1) then
         lingro = a
      else if (lambda .gt. 0) then
         xn = (2.0*w)**(1.0/3)
         call simp(func2,0.0,xn,int)
         aofxn = ((xn**3.0+2.0)**0.5)*(int/xn**1.5)
         if (a .eq. 1.64) write(*,*) xn,aofxn
         x = a*xn
         call simp(func2,0.0,x,int)
         aofx = ((x**3+2)**0.5)*(int/x**1.5)
         lingro = aofx/aofxn
         if (a .eq. 1.64) write(*,*) x,aofx,lingro
      else
         dn = func1(w)
         x = w*a
         lingro = func1(x)/dn
      endif
 
      end

      function func0(x)
      real x,func0
      func0 = 3/x + (3*((1+x)**0.5)/x**1.5)*log((1+x)**0.5-x**0.5)
      end

      function func1(x)
      real x,func1
      func1 = 1 + 3/x + (3*((1+x)**0.5)/x**1.5)*log((1+x)**0.5-x**0.5)
      end

      function func2(x)
      real x,func2
      func2 = (x/(x**3+2))**1.5
      end

      subroutine simp(func,a,b,s)
      implicit none
      real func,a,b,eps,ost,os,s,st
      integer jmax,j
      external func
      ost = -1.e-30
      os = -1.e-30
      eps = 1.e-5
      jmax = 20
      do j = 1,jmax
         call trapez(func,a,b,st,j)
         s = (4.*st - ost)/3.
         if (abs(s-os) .lt. eps*abs(os)) goto 3
         os = s
         ost = st
      enddo
      pause 'Too many steps.'
3     end

c
c     Trapezium rule from numerical recipes.
c

      subroutine trapez(func,a,b,s,n)
      implicit none
      real func,a,b,s,del,sum,x
      integer j,n,it,tnm
      external func
      save it

      if (n .eq. 1) then
         s = 0.5*(b-a)*(func(a)+func(b))
         it = 1
      else
         tnm = it
         del = (b-a)/tnm
         x = a + 0.5*del
         sum = 0.0
         do j = 1,it
            sum = sum + func(x)
            x = x + del
         enddo
         s = 0.5*(s+(b-a)*sum/tnm)
         it = 2*it
      endif
      end

      subroutine file_name_generate(filebase,ifile,filename)
      implicit none
      character*(*) filebase,filename
      integer ifile
      character*10 numeral
      integer len_trim,rindex

      write (numeral,'(i10)') ifile
      filename = filebase(1:len_trim(filebase))//'.'//
     &  numeral(1+scan(numeral,' ',.true.):len_trim(numeral))
      return
      end
      
      SUBROUTINE indexx(n,arr,indx)
      INTEGER n,indx(n),M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
      REAL a
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 13 j=l+1,ir
          indxt=indx(j) 
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END

